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ABSTRACT 

Motivated by the large number of extrasolar planetary systems that are near 
mean motion resonances, this paper explores a related type of dynamical be- 
havior known as "nodding". Here, the resonance angle of a planetary system 
executes libration (oscillatory motion) for several cycles, circulates for one or 
more cycles, and then enters once again into libration. This type of complicated 
dynamics can affect our interpretation of observed planetary systems that are 
in or near mean motion resonance. This work shows that planetary systems 
in (near) mean motion resonance can exhibit nodding behavior, and outlines the 
portion of parameter space where it occurs. This problem is addressed using both 
full numerical integrations of the planetary systems and via model equations ob- 
tained through expansions of the disturbing function. In the latter approach, 
we identify the relevant terms that allow for nodding. The two approaches are 
in agreement, and show that nodding often occurs when a small body is in an 
external mean motion resonance with a larger planet. As a result, the nodding 
phenomenon can be important for interpreting observations of transit timing 
variations, where the existence of smaller bodies is inferred through their effects 
on larger, observed transiting planets. For example, in actively nodding plane- 
tary systems, both the amplitude and frequency of the transit timing variations 
depend on the observational time window. 

Subject headings: planets and satellites: dynamical evolution and stability - 
planets and satellites: formation — planet-disk interactions 
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1. Introduction 

The current observational sample of extrasolar planets includes many systems with mul- 
tiple planets, and many systems have orbital period ratios that are close to integer values 
(e.g., Fabrycky et al. 2012). These systems are thus candidates for being in mean motion 
resonance (MMR), which represents a special dynamical state for a planetary system. In 
addition to the necessary period ratio, the other dynamical variables of a resonant system 
must allow one or more resonance angles (see below for their definitions) to execute oscil- 
latory behavior (e.g., Murray & Dermott 1999; hereafter MD99). One way to describe this 
requirement is that the resonant angle(s) must reside in a "bound state" within an "effective 
potential well" . Because of the special conditions required for a planetary system to reside 
in mean motion resonance, systems found in such states must have a constrained dynamical 
history. The relative fraction of planetary systems in mean motion resonance thus provides 
important information regarding planetary formation and early dynamical evolution. 

The dynamics of mean motion resonances is often more complex than indicated by 
standard textbook treatments. This paper explores one such complication called "nodding" , 
where the resonance angle librates for several cycles and then circulates for one or more 
cycles before returning to an oscillatory state. This paper addresses the problem using two 
complementary approaches, i.e., both full numerical integrations of the planetary systems 
and the construction of model equations obtained through expansions of the disturbing 
function. The expansion approach produces a large number of terms, and we identify the 
ones that allow for nodding behavior. These two approaches are in good agreement, and 
show that nodding often occurs when a small body is in mean motion resonance with a larger 
planet. The immediate goal of this paper is to obtain a better understanding of nodding 
in the context of three-body planetary systems. The over-arching goal is to provide a more 
detailed basis for interpreting observed systems that are found in or near resonance, including 
systems that exhibit transit timing variations (TTVs), which provide a means of detecting 
small bodies interacting with larger planets in transit (Agol et al. 2005). 

We note that previous dynamical studies, especially concerning asteroids in our So- 
lar System, have found behavior that is qualitatively similar to the nodding phenomenon 
explored herein. For example, asteroids near the 2:1 MMR with Jupiter are predicted to 
alternate between two modes of libration (Greenberg & Franklin 1975). In one mode, the 
longitude of perihelion for the asteroid librates about the longitude of perihelion for Jupiter; 
in the other mode, the aphelion of the asteroid librates about the longitude of conjunction. 
The asteroids are predicted to alternate between the two different modes. Similar apocentric 
resonances are found for the Hildas asteroid group (see Henrard et al. 1986; Ferraz-Mello 
1988). The phase space of these dynamical systems contain separatrices for both global and 
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internal (secondary) resonances (Morbidelli & Moons 1993); chaotic motion near the separa- 
trices can lead to crossing, and hence to alternating modes of libration (see also Michtchenko 
et al. 2008ab). In addition to appearing in planetary systems, this nodding behavior arises 
in other dynamical systems, including the driven, inverted pendulum (Acheson 1995). Ad- 
ditionally, nodding systems sometimes exhibit similar phase characteristics to other known 
dynamical systems, such as the the Duffing oscillator (Guckenheimer & Holmes 1983). 

Both nodding and transit timing variations can occur in planetary systems in or near 
mean motion resonance. Observed resonant and near-resonant systems provide important 
information about planetary systems. As one example, external perturbations can remove 
systems from resonance if the perturbations are large enough and if they act over a sufficiently 
long time (Adams et al. 2008; Rein & Papaloizou 2009; Lecoanet et al. 2009; Ketchum et 
al. 2011a); furthermore, such conditions can be realized in the circumstellar disks that 
form planets. As a result, planetary systems that are observed in resonance today must not 
have been greatly perturbed in the past, or they must have been subsequently influenced 
by significant dissipative interactions. On the other hand, resonances can also act as a 
protection mechanism, especially for close and massive planets. Indeed, the current sample 
of exoplanets includes systems that can only exist because of strong resonant interactions; the 
GJ876 system provides one such example (see Snellgrove et al. 2001; Beauge & Michtchenko 
2003; Kley et al. 2005). 

As another example, we note that entry into mean motion resonance is non-trivial: If 
the orbital elements of a planetary system are selected at random, the chance that the system 
resides in a mean motion resonance is relatively modest (even when the period ratio is chosen 
to be near the ratio of small integers). However, systems can evolve into resonance states 
through the process of convergent migration (e.g., Lee & Peale 2002), where, for instance, the 
outer planet migrates inward faster than the inner planet, and the two bodies subsequently 
move inward together. Even in this scenario, survival of the resonance can be compromised 
by overly rapid migration (Quillen 2006), and/or by turbulent forcing from the disk driving 
the migration (Lecoanet et al. 2009; Ketchum et al. 2011a). 

As one potential application of this work, nodding can affect our interpretation of TTVs 
(Agol et al. 2005). In this setting, unseen small bodies orbit outside observed transiting 
planets (usually Hot Jupiters). The smaller bodies affect the orbit of the inner, larger 
planets and lead to small variations in the timing of the transit events. This phenomenon is 
potentially a powerful method to detect (infer the presence of) smaller, otherwise unobserved, 
planets in such systems. Indeed, discoveries of this type have already been reported (e.g., 
Holman et al. 2010; Cochran et al. 2011), and many more are expected in the near future. 
However, the timing variations are largest when the smaller planets are in or near mean 
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motion resonance with the larger planet (e.g., Nesvorny & Morbidelli 2008), and such systems 
are susceptible to nodding as studied herein. Even without the complication posed by 
nodding, inferring the system properties from observations is a sensitive process (e.g., Veras 
et al. 2011). In any case, the results of this work will be useful for future interpretation of 
systems that exhibit transit timing variations. 

This paper is organized as follows. In Section [2j we study the nodding phenomena 
through numerical integrations of multiple planet systems that are near mean motion res- 
onance. This investigation shows that complex dynamical behavior, including nodding, is 
often present, and outlines the portion of parameter space where it occurs. For systems that 
exhibit nodding behavior, we then outline the corresponding effects on transit timing varia- 
tions. In Sectional we derive a class of model equations to describe nodding behavior. Here 
we expand the disturbing function for planetary interactions (e.g., MD99), keep the highest 
order terms, and identify the relevant terms that lead to nodding. The resulting model 
equations elucidate the dynamical ingredients required for nodding behavior to take place. 
We conclude, in Section 4, with a summary of our results, a discussion of their implications, 
and a brief description of future work. 



2. Numerical Study of Nodding 

2.1. Full 3-body Numerical Simulations 

This paper studies nodding of mean motion resonance angles for planet pairs that are 
near MMR. For simplicity, most of this work focuses on the 2:1 resonance. The term "nod- 
ding" here refers to a tendency to repeat a pattern of bounded libration for several cycles 
followed by one or more cycles of circulation. For intermediate times, the system exhibits 
behavior of MMR, but intermittent bouts of circulation may produce a cumulative net cir- 
culation of the resonance angle over many libration times. In the context of resonant angle 
phase trajectories, nodding can be described as motion near a separatrix in the phase space. 
The phase space for internal resonances contain one separatrix, whereas the phase space for 
external resonances can contain two distinct separatrices. The qualitative differences occur, 
in part, due to the existence of asymmetric external resonances, which arise when the orbital 
eccentricity for the test body becomes sufficiently large (Ferraz-Mello et al. 2003). The 
existence of asymmetric resonance strongly depends on the mass ratio (see Michtchenko et 
al 2008b); for example, for rri2 <C mi, all stable orbits are asymmetric (Beauge 1994). 

To carry out the study in this section, we numerically integrate the three-body gravita- 
tional forces using a Bulirsch-Stoer integration scheme. In addition to gravity, both general 
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relativistic corrections and stellar tidal damping are included in force calculations, but the 
inclusion of these additional forces (which are small) are not necessary to produce the inter- 
esting features explored in this work. Our system consists of a star with mass M* = 1M , 
a massive planet (m p = lMj up ), and a test body (m = 10 -6 Mj up ). We chose this particular 
test mass in order to minimize its influence on the planet's motion and to obtain the clearest 
dynamical signature of resonance angle nodding. However, nodding is also present for larger 
masses for the third body (in Section 2.4 we consider masses as large as m = 10M e ). The 
Jovian planet is placed in orbit with period T p — 1 year, and the test body is placed in 
orbit with initial period T — l/2orT = 2 years for studies involving internal and external 
resonances, respectively. We choose a benchmark test body eccentricity of e = 0.15, which 
is motivated by previous work (Ketchum et al. 2011a), and we choose from two values for 
the planet's initial orbital eccentricity, e p = 0.001 or e p — 0.1. 

To fully describe the initial configuration of the system, the initial orbital angles must 
be specified. We parameterize this study using the set of angles given by Azu - the relative 
alignment between the two orbits - and fo - the test body's true anomaly. Unless noted 
otherwise, all simulations begin with the planet and test body in conjunction. In an attempt 
to sample the available resonance angle phase space resulting from a choice of initial orbital 
elements {e p , e, a, Azu}, values of f spanning the full range — n to n in increments of n/100 
are sampled. Following this systematic approach, those states that have phase trajectories 
occurring near a separatrix, which ultimately lead to nodding, are easily found. 

The simulations are integrated for 10 4 years, a sufficient amount of time to capture 
several secular cycles for most initial configurations (for completeness, note that the param- 
eter space for external resonances contain small regions where the secular cycle's period is 
infinite - see Michtchenko et al. 2008b). The energy for a typical system that experiences 
no significant close encounters is conserved to better than one part in 10 10 . The planets' 
period ratio is monitored to confirm that the system remains nearly integer commensurate 
during the integration and hence that near resonance has not been compromised by a chance 
close encounter. The osculating elements for both bodies are recorded once per orbit of the 
inner planet. Since the libration timescales are of order ~ 30 — 100 orbits (albeit with large 
variations), this sampling cadence is frequent enough to resolve the behavior of the resonance 
angles. For 2:1 resonances, the resonance angles of interest are 

<f) = 2\ p -\-m, (1) 

for a test body internal to the planet's orbit, and 

= 2A - A p - w , (2) 

for a test body external to the planet's orbit, where A is the mean longitude and subscript 
p denotes the orbital elements belonging to the Jovian planet. And finally, the time rate 
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of change for the resonance angle, <f), is determined by quadratic interpolation, and used to 
construct resonance angle phase trajectories. 

2.2. Nodding Features for Near-Resonance 

The range of dynamical behaviors encountered in (or near) mean motion resonance is 
surprisingly rich (e.g., Michtchenko et al. 2008ab). Presented here is a small representative 
sample set of the simulations outlined above which display the main features found in these 
resonance states. 

We first consider systems where nodding does not occur (see Figured]). The figure fea- 
tures a system with a configuration and behavior deviating only slightly from the pendulum 
model of MD99. For this system, the perturbing planet is placed in orbit with semi-major 
axis a p = 1AU and eccentricity e p = 10~ 3 around a 1M star. A test particle of mass 
m = 10~ 6 Mj Up is set in a coplanar orbit with semi-major axis a = 0.63 AU and orbital ec- 
centricity e = 0.15, which places the two orbiting bodies near a 2:1 period ratio. The orbits 
are initially anti-aligned (Awq = tt) and the orbiting bodies placed in conjunction with the 
test particle near periapse - the system is prepared such that the planet's influence on the 
test body's motion is initially minimized, i.e., the two orbiting bodies cannot be separated 
further from one another during conjunction given this set of orbital elements. The top 
panel of Figure [1] shows the resonance angle, 0, from equation ([1]) in blue and the angle 
of apsides, Aw, in red. The resonance angle librates with a small amplitude of A<p ~ 0.1 
radians around the equilibrium = 0, while the apsidal angle circulates due to a prograde 
motion of the test particle's longitude of periastron - the planet's longitude of periastron, 
w p , does not move significantly. In the limit of small orbital eccentricities e and e p , the 
test body's true anomaly approximately coincides with the resonance angle at instances of 
conjunction. Thus, for this system, these dynamics depend mainly upon our choice for |/o| 
and are independent of our choice for Aw. The bottom panel of Figure [1] shows the phase 
trajectory for the resonance angle. This panel depicts a phase trajectory analogous to small 
oscillations of a simple pendulum, as expected. In this respect, the pendulum model of the 
circular restricted three body problem provides a sufficient model for resonances arising from 
orbital configurations of this type - although the planet's orbital eccentricity is non-zero, 
dynamical deviations from the pendulum model of the circular restricted three body problem 
incurred through small departures from circular symmetry are negligible. 

As |/o| increases from to tt, the amplitude of oscillations also increase until the con- 
junction line approaches apoapse, where the system reaches a separatrix and the resonance 
angle will circulate rather than oscillate. Figure |5] shows a system identical to that of Fig- 
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ure[T]in all aspects except that fo ~ 27r/3 radians, so that the resonance angle exhibits larger 
oscillations. The libration frequencies of the two systems however, are practically identical 
with frequency u ~ 0.03 yr _1 . 

For external resonance scenarios, the behavior is significantly different. Figure [3] shows 
the first example of nodding. In the figure, |/o| ~ n so initial conjunction is close to apastron 
of the test body's orbit, and the orbits are initially anti-aligned, Awq = tt. This set of initial 
orbital elements gives the maximum possible spatial separation between the planet and test 
body. The resonance angle 0, given by equation ([2]) and shown in blue in the top panel of 
the figure, seems to be attracted to one of two stable fixed points, with an unstable fixed 
point effectively located at <ft = tt. This is an example of an asymmetric resonance (Lee 
2004; compare with Callegari et al. 2004), and the existence of the two equilibrium points 
on either side of <ft — ft ar e due to a bifurcation occurring in the dynamics for sufficiently 
large test body orbital eccentricity (e.g., Michtchenko et al. 2008b). In this case, the 
system was configured near a different kind of separatrix than those encountered for internal 
resonances, and motion near this separatrix leads to nodding behavior that is unique to the 
external case. In general, the planar three-body problem has two degrees of freedom, with 
two proper frequencies: resonant and secular ones. As a result, oscillations of the problem 
(including nodding) must be combinations of the two proper frequencies. 

Note that there are several secular cycles shown in Figure [31 and during each secular cycle 
the resonance angle seems to choose one of three different libration modes - (a) resonance 
oscillations enclosing some point to the left of <ft = n, (b) oscillations enclosing some point 
to the right of = tt, or (c) oscillations enclosing all three points. As |/o| — » for similarly 
prepared systems, resonant oscillation amplitudes reach the maximum possible value A0 « tt 
radians before circulating. Encountered here is another separatrix, which is analogous to the 
separatrix for the internal resonance. In contrast to the internal resonance separatrix, the 
phase space trajectory for large amplitude oscillations here takes a different shape - the 
resonant angle's speed decreases (increases) as it approaches (departs from) the test body's 
apoapse location. Systems on such phase trajectories appear to nod once per one resonant 
libration, as depicted in Figure H] where fo ~ 7r/4. In this regard, the resonance angle moves 
as if it lives in a quartic potential, with a local maximum at <p = tt, two minima near <p = tt/2 
and (f) = 37r/2, and a maximum at cf) = 0. This feature of the motion near the outer separatrix 
of an asymmetric outer resonance can lead to a period increase or decrease by a factor of 
two for transit timing variations, a result presented later on in this paper. 

The nodding features become more prevalent as the perturber's eccentricity increases, 
where the resonance angle can begin accumulating a net circulation over longer times. Fig- 
ures [5] and [6] provide examples of the type of circulation behavior we observe in simulations. 
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These two figures also showcase the stark contrasts between nodding for internal and exter- 
nal resonances, respectively. In this comparison, the planet's orbital eccentricity is e p = 0.1 
- substantial enough to be well outside the circular restricted three body regime (where the 
pendulum model is strictly valid). Both examples show the tendency for Aw to circulate, 
on average, with a periodic component resulting from secular interactions, and, on average, 
the resonance angle circulates at the same rate. For times shorter than secular timescales, 
however, the resonance angle librates with some amplitude, not exceeding A0 = tt, around 
some equilibrium point. The point about which the resonance angle oscillates differs between 
the external and internal perturber cases. For the pendulum model of the circular restricted 
3-body problem, the equilibrium point is located at <fi = for internal resonances, while 
for external resonances it is = 7T. We stress, however, that the pendulum model is an 
over-simplification for the regime under consideration. 

Figure [5] shows the eccentricities (top panel), resonance angle (middle panel), and phase 
trajectory (bottom panel) for a system the same set of initial orbital elements as was used in 
Figure [U with the exception that here the planet's orbital eccentricity e p = 0.1. The system 
is prepared with anti-aligned orbits and the orbiting bodies in conjunction near the test 
body's apoapse. Recall from the above discussion that this particular configuration places 
the resonance angle near a separatrix in its phase space. The figure shows slightly more than 
3 complete secular cycles, each lasting w 3300 yrs, which corresponds to hundreds of libration 
times. During each secular cycle, the test body's orbital eccentricity (shown in black in the 
top panel) increases from its initial value of e = 0.15 up to e ~ 0.6 (where the two orbits 
intersect for a time), and then decreases close to its initial value. For times that are shorter 
than secular timescales and longer than libration timescales, the resonance angle undergoes 
large amplitude oscillation about the test body's periapse. In the previous case, when the 
planet's orbit was nearly circular, the motion of the test body's apsidal angle was steady 
circulation. In the present case where the Jovian planet's eccentricity is substantially larger, 
the azimuthal symmetry of the star-planet Keplarian system has been sufficiently broken 
and the evolution of the test body's osculating elements depends sensitively upon the orbital 
alignment (given by the apsidal angle, Aw). The red curve in the middle panel of Figure [5] 
shows the apsidal angle, which, over the course of one secular cycle, oscillates once very slowly 
about Aw = 7T, then very abruptly passes Aw = in the retrograde direction. Generally, as 
the apsidal angle approaches Aw = n (i.e., as the orbits approach anti-alignment), the test 
body's eccentricity decreases, and as the orbits rotate out of anti- alignment, the eccentricity 
increases. As the apsidal angle approaches Aw = 0, the orbits come into alignment, and with 
comparatively large orbital eccentricities, the planet exerts greater influence on the motion of 
the test body than in cases where it's orbit is nearly circular (see Batygin & Morbidelli 2011 
for a detailed analysis of secular dynamics). As a consequence, the resonance angle circulates 
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once or twice until the orbits rotate out of alignment, and the system enters into yet another 
secular cycle that deviates only slightly from the one just described. The test body can 
exhibit a wide range of eccentricity growth/decay during a secular cycle, and the duration 
of the secular cycle depends on the details of the initial 3-body configuration. However, the 
generic behavior of the osculating orbital elements described above during a secular cycle 
for any choice of Ati7 is robust. As long as the orbiting bodies are in conjunction with fo 
for the test body near apoapse, the resonance will reside near the separatrix in phase space 
and the system can experience nodding for states near mean motion resonance. 

For the external resonance shown in Figure |6l the system is initially configured with anti- 
aligned orbits Aw = n and conjunction occurring near the test body's apastron f ~ n. The 
test body's orbital eccentricity is sufficiently large to place the system near an asymmetric 
resonance, meaning there are generally two libration points. The exact locations of the 
libration points depend on the instantaneous orbital eccentricity (shown by the black curve 
in the top panel of the figure), but for ease of discussion, we take them to be ~ it ± 7r/2. 
During a secular cycle, the resonance angle <fi oscillates with an amplitude vr/4 < A0 < 7r/2. 
There are 6 full secular cycles shown here, each lasting ~ 1600 years. The test body's periapse 
circulates in the retrograde direction on average, so the apsidal angle Aro circulates. The 
nodding cycle, eccentricity growth/decay, and the apsidal angle circulation are all governed 
by the secular cycle and coordinate together to produce the specific nodding behavior of the 
resonance angle. As the two orbits' major axes rotate to become more perpendicular, the 
test body's eccentricity grows, and resonant oscillation amplitudes decrease, resulting in a 
tightening of libration. As the two orbits rotate so that their major axes become parallel 
(orbits either aligned or anti-aligned), the test body's orbital eccentricity decreases. The 
eccentricity attains a minimum value when the major axes are in alignment, reaching its 
smallest values when the orbits' periapses become aligned (as opposed to anti-aligned). As 
the test body's orbital eccentricity decreases, the libration points of the asymmetric resonance 
move closer to ir and the resonant angle's oscillation amplitudes generally increase, which in 
turn moves the phase trajectory toward one of the two separatrices. With the orbits near 
alignment (Aw ~ 0), the phase trajectory approaches the outer separatrix, and circulation 
becomes possible. When the orbits are nearly anti-aligned, the phase trajectory approaches 
the inner separatrix, which separates bounded libration around a single stable point from 
libration around both. Here, it is possible for the phase trajectory to jump across tc to the 
adjacent stable point. The number of times the resonant angle circulates or jumps back and 
forth between stable points during major axes alignment varies from secular cycle to cycle. 
However, the nodding cycles tend to mimic preceding cycles, and every so often abrupt 
changes will occur which again persist over multiple nodding cycles, until another abrupt 
change occurs. This behavior then repeats. 
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To summarize, we have studied both internal and external resonance scenarios using 
both nearly circular planet orbits (e p = 10~ 3 ) and an orbit that departs from circular (e p = 
0.1). We took a benchmark value for the test body's initial orbital eccentricity of e = 0.15, 
but this value varied with the secular cycle, reaching values as high as e ~ 0.6 and as 
low as e ~ 0.01. Variations to the test body's orbital eccentricity increased in simulations 
which included larger planet eccentricity. On secular timescales, the 2:1 resonance angle 
may exhibit bouts of circulation with regular libration in between circulation events (the 
phenomenon we call nodding). Nodding is common when the planet's orbital eccentricity is 
sufficiently large, where we find an approximate threshold of ~ e p > a few xl0~ 2 . Nodding 
systems with larger orbital eccentricities are able to obtain greater net circulation than 
systems where both orbital eccentricities are small. However, the relative orbital alignment 
and the test body's true anomaly at the moment of conjunction affect both the period of 
the nodding cycle and the libration width during intermediate times. Both internal and 
external resonances exhibit nodding, but prominent qualitative differences in the nodding 
signatures between the two configurations exist. Because of our choice of initial test body 
orbital eccentricity, asymmetric resonances are found for external cases, where libration 
occurs around points = (l±5/2)n with 5 1 for e ps 0.15 (e.g., Lee 2004; cf. Callegari et 
al. 2004). Consequently, for external resonances, libration amplitudes are typically smaller 
(A0 < 7r/4) for large e p . For nearly circular planet orbits, large amplitudes (A<f> < n) may 
persist for external resonances, and appear to librate about = ir. However, the resonance 
angle's motion slows as it passes = ir, which is distinctively different from large amplitude 
oscillations for the internal resonance case. 



2.3. Dynamical Map of Nodding in the Hot Jupiter Problem 

In this section we explore the nodding phenomenon for Hot Jupiters. Because of the 
existence of the asymmetric resonance and the additional separatrix that comes along with it, 
the external resonance is particularly interesting in the context of the Hot Jupiter problem 
(where smaller bodies could be found in outer orbits - see Ketchum et al. 2011b). For 
completeness, we explore how the resonance angles of the external 2:1 near resonance are 
affected by the external test body's orbital eccentricity and by the orbital alignment between 
the planet and the test body. To perform this study, we use typical parameters of Hot 
Jupiter systems, where a Jovian planet (with mass m p = lMj up , orbital eccentricity e p = 
0.04, and orbital period T p = 4 days) orbits a 1M star with a small planet (with mass 
m = 1M© and orbital period T = 8 days) orbiting external to the Star- Jovian system. We 
take orbital eccentricity values for the test body between —3 < log 10 e < —1/3 in equally 
spaced logarithmic increments. For each value of the test body's eccentricity, we perform an 
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Fig. 1. — Internal 2:1 near-resonance with periaspses initially anti-aligned (Aw = it), test 
body near periastron at conjunction (/ ~ 0), and e p = 1CT 3 . All angles are in radians and 
time is in years. The top panel shows the resonance angle, <ft (blue curve), and the angle 
of apsides, Azu (red curve). In this example, resonance angle nodding does not occur. The 
bottom panel shows the phase trajectory of the resonance angle during the full 10 4 years 
of simulation. This system is close to the text book example of mean motion resonance 
described by the pendulum model (see Murray & Dermott 1999). The libration amplitude 
is small (A0 ~ 0.1 rad) and the equilibrium point for oscillations is at = 0. 
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Fig. 2. — Internal 2:1 near-resonance with periapses initially anti-aligned (Aw = 7r), test 
body closer to apastron than to periastron during initial conjunction with a true anomaly 
f m 2tc/3, and the perturbing planet's orbit nearly circular with an eccentricity of e p = 10~ 3 . 
All angles are in radians and time is in years. The top panel depicts the resonance angle 
(blue curve) and the angle of apsides (red curve). The resonance angle undergoes large 
amplitude oscillations of A(f> > 37r/4 due to the test body's large angular displacement 
from periastron of its orbit during the initial conjunction with the perturbing planet. In 
this example, resonance angle nodding does not occur. The bottom panel shows the phase 
trajectory for the resonance angle over the full 10 4 years of the simulation, which oscillates 
about (modulo 2tt) and doesn't circulate at anytime. 
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Fig. 3. — External 2:1 near-resonance with periapses initially anti-aligned (Aw « tt), test 
body near apastron at conjunction (/o ~ tt), and a nearly circular planetary orbit with 
eccentricity e p = 10~ 3 . All angles are in radians and time is in years. The top panel 
shows the time series of the test body's eccentricity (black curve - taking values between 
0.14 < e < 0.18) and the planet's eccentricity (blue curve - nearly constant at e p ~ 10~ 3 ). 
The middle panel shows the resonance angle, (blue curve), and the angle of apsides, Aw 
(red curve). The resonance angle jumps between two points of attraction located near 7t/2 
and 3vt/2 (mod 2ir) - sometimes librating around one of these two points many times (e.g., 
for times between 1000 < t < 3000) before jumping to the opposite point of attraction - this 
is one form of nodding. The bottom panel shows the resonance angle's phase trajectory for 
the entire 10 4 years of simulation. This trajectory traces out a complex dance around the 
inner separatrix of the phase space of the external 2:1 asymmetric resonance. 
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Fig. 4. — External 2:1 near-resonance with periapses initially anti-aligned (Aw ~ it), test 
body displaced ~ 40° from periastron at conjunction (/ ~ tt/4), and a nearly circular 
planetary orbit with eccentricity e p = 10~ 3 . All angles are in radians and time is in years. 
The top panel shows the resonance angle, (j> (blue curve), and the angle of apsides, Aw (red 
curve). The resonance angle repeatedly jumps between two points of attraction located near 
7r/2 and 3tt/2 (mod 2tt) - this behavior is one form of nodding. The bottom panel shows 
the resonance angle's phase trajectory for the entire 10 4 years of simulation. The resonance 
angle's motion temporarily slows as it passes apoastron of the test body's orbit (f> « — ir, 
and accelerates producing sharper peaks in the resonance angle signal of the top panel, in 
contrast to the internal resonance case (compare with Figure EJ. 
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Fig. 5. — Internal 2:1 near-resonance, periapses initially aligned (Aw w 0), test body near 
apastron during initial conjunction (/ ~ tt), and the external planet's orbit with eccentricity 
is e p = 0.1. All angles are in radians and time is in years. The top panel shows the eccentricity 
time series for both the test particle (black curve) and the planet (blue curve). The test 
body's eccentricity reaches values in excess of e > 0.6 and as low as e < 0.1, the planet's 
orbital eccentricity does not vary significantly. The middle panel shows the time series for the 
resonance angle (in blue) and Aw (in red). The resonance angle librates about 0^0 with 
a varying amplitude tt/2 < A<p < tt for large stretches of time, and suddenly circulates once 
(at time close to t ~ 6500 years) or twice (at time t pa 3000 years). This is a form of nodding. 
The moments in time where the nodding events occur are correlated with the times that the 
test body's orbital eccentricity is comparable to the planet's orbital eccentricity, e ~ e p , and 
to times when the orbits are aligned (the angle of apsides Aw ~ modulo 2tt). The bottom 
panel shows the corresponding phase trajectory of the resonance angle, (ft vs (p. 
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Fig. 6. — External 2:1 near-resonance, periapses anti-aligned (Aw « it), test body near 
apoastron during initial conjunction (/ ~ n), and planet on an eccentric orbit with e p = 0.1. 
All angles are in radians and time is in years. The top panel shows the eccentricity time 
series for both the test particle (black curve) and the perturber (blue curve). The planet's 
eccentricity remains nearly constant during the simulation, while the test body's eccentricity 
reaches values up to e ~ 0.4 and down to e ~ 0.05. The middle panel shows the time series 
for the resonance angle (in blue) and Aw (in red). The bottom panel shows the phase 
trajectories of the resonance angle, vs <j). 
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ensemble of 200 similarly prepared systems, each with a slightly different initial apsidal angle 
between —tt< Aw < tt in equally spaced increments. Both the Hot Jupiter and the smaller 
rocky planet begin each simulation located at the periastron of their respective orbits (not 
necessarily in conjunction with one another). The parameter space outlined by the above 
conditions contains 28,560 points, each of which are integrated for up to 200,000 orbits of 
the Jovian planet (~ 2000 years). During each individual simulation, we keep track of the 
displacement of both the resonance angle and the angle of apsides Aw, as well as the total 
number of times the resonance angle passes the test body's periastron (the latter provides a 
rough measure of the number of circulation events that take place during the simulation). 

The results of the survey are shown as a dynamical map in Figure Each pixel 
in the figure is composed of an admixture of the three colors red, green, and blue. The 
amount of each color within a given pixel represents different behavioral characteristics of 
the angles of interest: (i) Red measures the resonances angle's final angular displacement, (ii) 
Green measures the angle of apsides final angular displacement, and (Hi) Blue measures the 
resonance angle's final angular displacement in comparison to the total number of circulation 
occurrences during the given simulation. We use the ensemble averages for each category to 
gauge the overall shade of each color of the survey. 

Admittedly, our metric for defining each color's particular shade in each pixel is some- 
what arbitrary, so we refrain from providing specific details about it here. However, in 
general our metric produces a bright color for smaller deviations from the initial state in 
comparison to the ensemble's average deviation, and a darker shade corresponds to a larger 
deviation than average. We provide a brief interpretation of the prevalent color combinations 
presented in the figure as follows: RED: strictly librates (no circulation events) and Aro 
circulates; YELLOW: Both and Aw librate (minimal circulation); GREEN: Aw librates, 
but circulates quickly; PURPLE: circulation events encountered, but direction of circu- 
lation is erratic, and Aw ciruclates; CYAN: relatively few circulation events for both and 
Aw and direction of circulation is consistent; BLACK: both and Aw circulate rapidly; 
BLANK: either (i) premature termination of the simulation due to a scattering or collision 
event for the test body, or (ii) a period ratio between the planet and test body that deviates 
by more than 10% once the simulation reaches its time limit. 

Although the details appearing in the Figure [7J are intricate and rich, there two main 
features that we want to emphasize here. [1] The existence of the large, red islands appearing 
near (log(e) , Aw ) = (— 0.5,7r/2) and the absence of the red islands for values log(e) < 
— 1 is in agreement with previous findings that for m 2 <C m 1; stable oscillations orbits 
are asymmetric (Beauge 1994). These large red islands are regions of parameter space 
where the resonance angle does not circulate at any point in time during the 200,000 orbit 
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simulation. [2] The diagram shows an apparent pitchfork bifurcation occupying this mapping. 
For low test body eccentricities, rapid circulation occurs for initial apsidal angles, indicated 
by the thin horizontal tracks of dark (almost black) pixels located around Aw = ±7r/2. 
(For completeness, note that the secular angle is not defined in the limit e — > 0; here the 
eccentricity axis is presented on a logarithmic scale so that e > 0.001.) As the eccentricity 
is increased, these tracks converge to a single horizontal line at the center of the map, where 
at e w 0.1 rapid circulation of both and Aw occurs for initial apsidal angle Azu 0. 
This trend in the dynamics of the resonance angles, as the exterior test body's eccentricity 
is increased, is consistent with past studies (e.g., Lee 2004) and, perhaps, reveals further 
details about the onset of the asymmetric exterior 2:1 resonance. 
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Fig. 7. — Dynamical map of the resonance and apsidal angles' relative amounts of accumu- 
lated circulation during simulations lasting 200,000 orbits (~ 2000 years) for the external 
2:1 near-resonance. All angles are given in radians. Each pixel in this figure represents 
one simulation with a given initial apsidal angle, Aw (vertical axis), and test body orbital 
eccentricity, log 10 e (horizontal axis). This figure is a mosaic of 28,560 different simulations 
performed as described in the text. Each pixel is colored with an admixture of red, green, and 
blue. The particular shade of each color signifies the relative amount of circulation, hence 
nodding, using three different metrics - see text for full description. This set of simulations 
was performed under the context of a Hot Jupiter scenario, i.e., the planet has eccentric- 
ity e p = 0.04, mass m p = 1 Mj Up , and period T p = 4 days, and the test body has mass 
m ~ 1M© and period T = 8 days, both orbiting a 1M Q star. Each simulation begins with 
the planet and test body located at periastron of their respective orbits. This figure shows 
what looks like a pitchfork bifurcation in this parameter space (perhaps due to the transition 
from symmetric to asymmetric resonances following from left to right across map). For low 
test body eccentricities, rapid circulation occurs for initial apsidal angles Aw = ±7r/2, and 
as the initial eccentricity increases, these tracks converge in the center, where, at e ~ 0.1, 
rapid circulation occurs for states with an initial apsidal angle Aw ~ 0. 



-20 - 



2.4. Transit Timing Variations in the Presence of Nodding 

One potential channel for discovering terrestrial sized exoplanets involves observing 
variations in the transit times for a Hot Jupiter (e.g., Agol et al. 2005). These transit timing 
variations (TTVs) will be largest if the body responsible for the variations is in (or near) 
a mean motion resonance with the Hot Jupiter. Of the multiple-planet systems known to 
date, many of the adjacent planet pairs have period ratios that are close to small integer 
ratios, with 2:1, 3:2, and 3:1 being the most populated, in order of decreasing frequency (e.g., 
Fabrycky et al. 2012). These observations suggest that a planet, whose presence produces 
transit timing variations of a transiting Hot Jupiter, is more likely to be near MMR, not 
deep in MMR with the Hot Jupiter. On the other hand, planet-pairs could still be deep in 
resonance while the ratio of their orbital periods deviates from exact resonance by as much as 
20% (e.g., Batygin & Morbidelli 2012). As shown in the previous section, planetary systems 
near MMR can exhibit distinctive and interesting behavior — nodding of the resonance 
angle. This phenomenon can lead to interesting signatures in TTVs, which are considered 
in this section. 

As outlined above, in order for nodding to arise, the larger planet must have nonzero 
eccentricity. However, Hot Jupiters are subject to tidal forces, which act to circularize their 
orbits. On the other hand, the observed orbits of Hot Jupiters are often eccentric: The 
current sample includes 242 planets with semimajor axes a < 0.1 AU (Exoplanets Data 
Explorer 2012) and 79 of these planets have orbital eccentricity in the range 0.02 < e < 0.6 
(and 41 orbits have e > 0.10). Some fraction of the Hot Jupiters are thus expected to 
undergo nodding (if smaller companion planets are also present). 

The Hot Jupiters could have nonzero eccentricity for several reasons: [1] The tidal 
circularization time for 4-day orbits is typically r C j rc ~ 1 Gyr (Goldreich & Soter 1963; Hut 
1981), and observed systems have ages of 2 - 5 Gyr. As a result, the initial eccentricity is 
decreased by a factor of "only" ~ 7 — 150. If, for example, the starting eccentricity e = 0.30, 
the observed eccentricity e 0.002 — 0.04; younger systems with higher starting eccentricities 
can have even larger values. [2] The tidal circularization time grows with a high power of 
semimajor axis, r circ ~ a 13//2 (Hut 1981), so that planets with slightly larger orbits have 
much longer circularization times and hence larger eccentricities. [3] The time required for 
tidal circularization depends on the tidal quality factor Q (Goldreich & Soter 1963), which 
depends on the internal structure of the planet. Since observed planets are inferred to have 
a wide range of structures for a given mass (e.g., Bodenheimer et al. 2003), we expect a 
wide range of Q values. [4] In multiple planet systems, massive distant planets can drive 
Hot Jupiters to have nonzero eccentricities (Adams & Laughlin 2006). 

To demonstrate the possible effects of nodding on transit timing variations of a Hot 



-21 - 



Jupiter, we perform additional simulations. We place two planets in orbit around a 1M & 
star - a lMj up planet in a 4 day orbit and a 1M® planet in an 8 day orbit. The periastra are 
approximately anti-aligned Ace « 7T, and the planets are initially placed near conjunction 
with the Hot Jupiter at periastron. The orbits are set to be coplanar. The initial eccentricity 
of the smaller planet is taken to be e = 0.15, whereas the eccentricity of the Hot Jupiter 
has values e p = 0.1 and e p = 0.01. For each choice of e p , we then performed numerical 
integrations using the integration scheme described above and recorded the times that the 
center-to-center displacement between the star and Hot Jupiter were parallel to some arbi- 
trary line of sight. Neither light travel time nor transit duration variations are corrected for 
during the simulations - a more rigorous study should take into account the light travel time 
from the star to the Hot Jupiter, and should consider the surface-to-surface grazing chord 
between the star and Hot Jupiter as the line-of-sight rather than the center-to-center line 
(Veras et al. 2011). However, these are higher order corrections and are not considered here. 

Figure |8] shows the results of one such simulation for the larger eccentricity e p = 0.10 of 
the Hot Jupiter. The resonance angle, shown in blue on the top panel, begins in one mode 
of oscillation, is caught in libration while being passed back and forth between fixed points 
near —it/2 and — 3ir/2, and then settles on the — 3n/2 fixed point for nearly 20 librations. 
At t — 10 years (~ 1000 orbits), the periastra have rotated into near alignment Aw w 0, and 
the resonance angle undergoes a different type of nodding than exhibited at the beginning, a 
mode which is characterized by free circulation across the unstable point located at = — 2n. 
Near the t — 12 year mark, the resonance angle enters into a mode of libration around a single 
fixed libration point at ~ — 57r/2. This libration undergoes around 20 oscillations, at which 
time the cycle starts over from the beginning and is repeated several times. The particular 
cycle illustrated here, which produces the nodding pattern, is similar to that exhibited by 
the system represented in Figure [6j the initial orbital configurations are nearly identical. 
However, the mass of the smaller body for this case, with m = 1M®, is much larger and the 
two planets are located much deeper inside the stellar potential than in the previous case, 
deep enough that stellar damping plays a role in the dynamics. These differences contribute 
to producing the different signatures of the nodding patterns for the two cases. 

The transit timing variations described above and shown in Figure [H] exhibit many 
different modes of oscillation. In addition, this figure shows that it is plausible that the 
transit timing variations closely follow the time derivative of the resonance angle. When 
the resonance angle is nodding between the fixed points on either side of 7r, the derivative 
decreases as passes it, and then increases before reaches the stable point opposite n and 
reverses direction. This type of motion qualitatively matches the transit timing variations 
calculated during the same period of time. When the resonance angle increases, the transit 
timing variations are positive; when the angle decreases, the TTVs are negative; and when 
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the angle passes through ir, the angle slows down and the TTVs decrease. This complicated 
behavior is depicted by the 'double-peaked' patterns in Figure [8] between t = — 2 years and 
t = 20 - 25 years. 

The bottom panel in Figure [8] shows the Fourier Transforms (FT) of the TTV signal in 3 
year windows which encapsulate times corresponding to the three different modes of libration 
exhibited by the resonance angle's nodding cycle. The three different power spectra in the 
bottom panel are color coded and correspond to the regions in the TTV signal highlighted 
with the same color. Accordingly, the dashed blue curve is the FT for the portion of the 
TTV signal in the middle panel that is highlighted in blue, near t ~ 5 — 8 years. During this 
window, the resonance angle librates around a single fixed point, and the peak frequency 
given by the FT is / ~ 3 yr -1 . The red FT curve is from the TTV signal in the time 
window of t ~ 10 — 13 years, which corresponds to a dynamical state of the system which 
exhibits a net resonance angle circulation with aligned orbits. The peak frequency of the 
red FT curve is / ~ 1.4 yr -1 . Finally, the green FT curve from the TTV signal between 
t ~ 21 — 24 years has a broad peak near frequencies / ~ 0.5 — 0.8 yr -1 . Thus, for this 
particular system, observations of TTVs over a three year time interval can yield a range 
of fundamental frequencies that span a factor of ~ 4; the result depends on the three year 
window for which the transit times are observed. 

Figure |9] shows the time evolution for a planetary system that is identical to that shown 
in Figure El except that the Hot Jupiter is taken to have smaller eccentricity e p = 0.01. The 
behavior of the resonance angle is much the same as before, i.e., the transit timing variations 
calculated for the Hot Jupiter in this case contain all of the same qualitative features dis- 
cussed in the previous example. For completeness, we have run additional scenarios where 
the mass of the smaller planet takes on different values. The resulting time evolution of the 
resonance angles (not shown) are, again, qualitatively similar. As the mass of the smaller 
body varies, however, the amplitude of the transit timing variations change in proportion, 
as outlined below (see also Agol et al. 2005; Veras et al. 2011). 

Before leaving this section, we derive an approximate relation between TTV amplitudes 
(for the Hot Jupiter scenarios outlined above) and the time rate of change for the resonance 
angle <fi. The applicable resonance angle for this situation is given in equation (T5]). The 
second time derivative for the resonance angle is given by 

4>= (P + Q)(n + '£) S) -p(n + e) j -qiu (B 
~ (p + q)n e - phj . 

Lagrange's planetary equations of motion provide us with the time rate of change for the 
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mean motion for both the Hot Jupiter and the Super Earth in this problem, 

3 

n® = — r (p + q)d<j>Ki , (4) 
a e 

3 

= —pd^R, E , (5) 

a j 

where TZj and TZe are the disturbing functions for an internal and external perturber, re- 
spectively. Using a simple, time averaged form for the internal disturbing function, 

TYl ' 

IZi = [aS + -Fie® cos 0] , (6) 

and the external disturbing function, 



a 3 



a® 

where S is the secular contribution to the disturbing function, we substitute these forms into 
equations (HJ and (111) and get a relation between and hj, 

P + q J 7 ! 

where /i = m^/rrij is the ratio of planet masses. In the case of a 2 : 1 resonance, a ~ 2~ 2 / 3 , 
and the portion of equation ([8]) that depends on a simplifies to 

3 /31 — 2a 2 

a -=- — a — ; — — - = a . (9) 

7i a/ 31 - (2a)- 1 1 ; 

Substituting this simplification into equation and combining the result with the approx- 
imate form for the second time derivative in equation ([3]), we find an expression relating the 
time rate of change for the Hot Jupiter's orbital period to the second time derivative of the 
resonance angle, 

^ « *rn^ > ( 10 ) 

where we have made use of the relation h = —2ttT 2 T in this last line. Neglecting time 
variations in a, equation ffTO]) can be directly integrated (over one orbit of the Hot Jupiter) 
to give us the desired relation between transit timing variations, AT, and the time derivative 
of the resonance angle, <f), 

liTf . 



To demonstrate the accuracy of equation (fTTj) . we take Figure [8] as an example. At peak 
TTV in Figure [HI (near the 5 year mark in the TTV panel of the figure), the resonance angle 
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is changing at a rate of ~ 12 rad yr -1 ~ 3.8 x 10~ 7 rad sec -1 . The Hot Jupiter's period is 
roughly 4 days (~ 3.5 x 10 5 sec), and the outer planet has mass m = 1 M e (/i w 0.003). 
Substituting these quantities into the right hand side of equation gives a timing variation 
AT « 14 sec, a result which is in good agreement with the figure. This result shows that the 
TTV signal can reveal intricate details about near-resonance nodding if such a system can 
be found in nature. The above derivation also confirms that the strength of the TTV signal 
increases with the mass of the exterior planet (Agol et al. 2005). We have run additional 
simulations with varying masses for the outer planet, up to m = 10 M e , and find at the 
amplitudes of the TTVs vary in proportion to the mass (as expected). In addition, equation 
ffTTj) shows that the TTV signal increases as the square of the Hot Jupiter's orbital period 
and one power of <fi. Since itself is inversely proportional to the period, the TTV amplitudes 
are directly proportional to the orbital period. 
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Fig. 8. — Results for a Hot Jupiter in a 4-day orbit, with eccentricity e = 0.10, and in an 
external 2:1 near-MMR with an Earth-mass planet (m m 1M®), both orbiting a 1M star. 
The top panel shows the time series for the resonance angle <fi (blue curve) and the angle 
of apsides Aw (red curve) for the first 25 years (~ 2000 Hot Jupiter orbits). The initial 
configuration of this system was chosen to obtain nodding behavior of the resonance angle. 
The middle panel shows transit timing variations (in seconds) for the Hot Jupiter during 
the simulation. Three separate 3-year windows highlight the TTV data (corresponding to 
different nodding modes of the resonance angle) in the colors (from left to right) blue, red, 
and green. Fourier Transforms for each highlighted portion are shown in the bottom panel, 
and are color coded to match their corresponding progenitive TTV signal. The bottom panel 
demonstrates that the fundamental frequency of the TTV signal for a system undergoing 
nodding motion can vary by a factor ~ 2 (either way) depending on resonance angle nodding 
effects. 
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Fig. 9. — Results for a Hot Jupiter in a 4-day orbit, with eccentricity e = 0.01, and in an 
external 2:1 near-MMR with an Earth-mass planet (m ~ 1M®), both orbiting a 1M star. 
The top panel shows the time series for the resonance angle (f> (blue curve) and the angle of 
apsides Ave (red curve) for w 35 years (~ 3000 Hot Jupiter orbits). The initial configuration 
of this system was chosen to obtain nodding behavior of the resonance angle. The middle 
panel shows transit timing variations (in seconds) for the Hot Jupiter during the simulation. 
Three separate 3-year windows highlight the TTV data (corresponding to different nodding 
modes of the resonance angle) in the colors (from left to right) blue, red, and green. Fourier 
Transforms for each highlighted portion are shown in the bottom panel, and are color coded 
to match their corresponding progenitive TTV signal. The bottom panel demonstrates that 
the fundamental frequency of the TTV signal for a system undergoing nodding motion can 
vary by a factor ~ 2 depending on resonance angle nodding effects. 
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3. Derivation of Model Equations for Nodding 

In this section, we turn our attention to a Lagrangian formalism for the scenarios out- 
lined in section 2. We analyze the equation of motion for the resonance angle in the restricted 
3-body problem in order to identify the terms that encapsulate the dynamics of nodding. 
We choose a Lagrangian formulation instead of a Hamiltonian treatment purely for ease of 
physical interpretation. Consider three point masses that follow a hierarchical arrangement, 
m t <C m p <C m*, where we denote orbital elements belonging to the test body with subscript 
t and those belonging to the planet with p. The three bodies are spatially separated so that 
m t and m p are gravitationally bound to m*, but not bound to one another. Furthermore, 
suppose the orbits for m t and m p are coplanar, so they orbit m* in the same plane. In 
what follows, first we provide a general treatment of the equations where we do not specify 
whether the planet's orbit is interior or exterior to the test body's orbit, but will later focus 
on the internal perturber case where a p < a t . The system sketched in the above outline is the 
planar restricted 3-body problem. This coming derivation follows a mathematical approach 
similar to that used in the text by Murray & Dermott (1999). 

The argument angles (and hence the resonance angle), 0, for the disturbing function in 
the planar 3-body problem are given by a linear combination of the Eulerian angles 

= 3 At + jpAp + k t w t + k p Wp , (12) 

where j t , j p , k t , and k p are integers, which sum to 0, j t + j p + k t + k p = 0. For such a system 
in resonance, \j p \ : \j t \ is a ratio of small integers, and j p and j t are of opposite sign. From 
simulations, we found that the angle that tends to librate first is the angle corresponding 
to that of the circular restricted 3-body problem, i.e., where k p = in the case where the 
perturbing planet's orbit is strictly circular. This is because the strength functions associated 
with the resonance angle (besides those of the secular angles) are proportional to k p factors 
of the planet's orbital eccentricity e p and k t factors of the test body's orbital eccentricity e t . 
The perturbing planet and the primary mass (star) form a nearly Keplarian system. As a 
result, we assume that X p « n p in this derivation. From this point forward, we note that 
orbital elements appearing without a subscript belong to the test body. 

We find the first time derivative for the resonance angle 

<t> = jt (n + e) + j p n p - (j p + j t )w , (13) 

where the quantity 

e = e + ni (14) 

is introduced to avoid explicit time occurrences in the equation of motion (Brouwer & 
Clemence 1961) and where h in equation ( I14p may originate from both damping effects 



-28 - 



and the disturbing function. In the above derivation, it is important to note a few things. 
From this point forward we use the following notation; partial derivatives written as d q op- 
erate only on explicit occurrences of the variable q, and otherwise vanish. Using a lowest 
eccentricity order form of Lagrange's planetary equations of motion for e and w, we can 
write 

w = -^rd e K , (15) 
na z e 

and 

2 e 2a e 2 

i = — d a n + 7-^d e n = s A —d a n + -w , (16) 

na 2na z na z I 

where TZ is the disturbing function (for either an interior or an exterior perturber), and where 
we introduce Sa with Sa = {1, — 1} for an internal or external perturber, respectively; this 
sign function arises from a change of variables from at to a, i.e., ad a = —SA&d a . We can 
rewrite equation (|T3l) as 

4> = (j t + na SB V{K)) n + j p n p , (17) 

where fi = m p /m* is the mass ratio between the planet and the star, we introduce the 
quantity Sb where Sb = {0, 2} for an internal or external perturber, respectively; we have 
also introduced the linear differential operator V defined as 



V = A 1 (h a d a - h e d e) 



In equation ( TL8l) . A = Qm p /a p is the pre-factor to the disturbing function containing dimen- 
sionful constants, and we introduce the quantities h a , h e , and defined by 



h __ (jt+jp)P 

e ~ ae ' (19) 
jt c 2 

and j3 = 1 . 

jt + 3 P 2 

Note that when substituting Lagrange's planetary equation of motion for e in equation ( TL6l) . 
we have kept contributions from variations to the disturbing function from a, whereas some 
treatments will ignore these contributions (e.g., MD99), equivalently asserting that d a ~ 0. 
However, this contribution can be substantially larger (nearly 40 times greater) than the 
contribution to e from the we 2 /2 term. The form in equation (JTTJ) makes it easier to take 
the second time derivative en route to the resonance angle's equation of motion, 



[jt + /ia SB V(n)] h + fina SB S B ^V{TZ) + ^ {V{TZ)} 



(20) 
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Using the form of the operator V given in equation f)T8]) . we can apply the chain rule for the 
time derivative in the last term of equation f )20|) to obtain 

d 



dt 



[V(TZ)\ =V(K)+V(K) . 



(21) 



The time derivative of the V operator, for constant primary and secondary masses, and 
constant semi-major axis for the secondary mass, is given by 



V{K) = A- L h 



J 



and 



2n 



dJZ , with (3* 



V(jl) = eV{d e K) + S A ~aV{d a Tl) + wV(d m TL) + <pV{d^Tl) 

3 n 



(22) 



(23) 



Combining the previous two results with equation (|20|) gives 



n 
n 



J P n p + S A -fina SB {S B V(K) + aV(d a TZ) + k^h e d e K) 



+ fxna 



6 - (eV(d e 1l) + ^ l h e j d e TZ ) + <j>V{d+K) + wVid^ll) 



(24) 



Equation (|24|) is arranged so that disk damping effects can be easily included into the equa- 
tion of motion - all time derivative operations of orbital elements have been collected and 
accounted for and occur as first order time derivatives. With this form, one can simply 
substitute Lagrange's planetary equations of motion where necessary and remove all implicit 
instances of time from the equation of motion, i.e., all remaining operators act only as partial 
derivatives with respect to the orbital elements, independent of time. Along these lines, one 
could modify the time derivatives of the orbital elements to include external forces (e.g., 
damping) 



(25) 



n = (n)dist + Wdamp J 

e = (e)dist + (e)damp , 

where (g)dist represent variations to the orbital element q due to the disturbing function, 
and (g)damp are variations caused by external forces due to interactions with, for instance, a 
circumstellar disk. 



3.1. Specialization to the Case of External resonances 

We will now focus an external p + q : p resonance. For a given angle (f) from equation 
( TT2"]) . The simplest time averaged disturbing function (to second order in eccentricity) for an 
internal perturber takes the form 

IZ^ = A [a (Si + S 2 cos (w — zo p )) +1Zd cos 4>] , (26) 
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where 



A 



5 1 = (e 2 + e 2 p ) f sA (a) , (27) 

5 2 = ee p f St2 (a) , 

TZ D = e q J r i(a) where 7i(a) = (af d (a) + /»(a)) , 

(MD99). Under this scenario, the following quantities will have values j t = (p + q), j P = —p, 
Sa — 1, and 5^ = 0. Using the damping forms given above in equation fl25|) . we rearrange 
equation (|2"%j) algebraically and use the time averaged form of the disturbing function in 
equation f l26l) . and parse the equation of motion into generalized reoccurring operators, 
trigonometric functions of <p and Aw, and damping terms. Note that, during this expansion, 
we expand the equation for e to the second lowest order in eccentricity to retrieve a form 
that utilizes the parameter (3. This procedure results in an equation of motion of the form 



— sin i 



Afpn p + (Af + fmAV(TZ D ))<j) + J(o&l) + J(aS 2 ) cos (Aw) + J(K D ) cos , 



— sin (Aw) [fj,nwAV(aS 2 ) + K,(aS\) + K,(aS 2 ) cos (Aw) + JC(R,£>) cos< 
+ Af(pn p + 4>) + J(aS x ) + J(aS 2 ) cos (Aw) + J(H D ) cos* 

Here, new operators were introduced to compactify notation 

J(Q)=NV n (Q)+8 4) V e (Q) , 
K(Q) = S^V e (Q) , 
J(Q)=MV n (Q)+SV e (Q) , 

with 

V n (Q) = ^fin [aAV(d a Q) + Kd e Q] , 

P* 



(28) 



(29) 



V e (Q) = fin 

where Q is some function of eccentricities and a given in equation (1271) . 



eAP(9 e Q) + /i e ^-9 e Q 



(30) 



N= -3 C(p + q)K D 



£-vj = —C—=S 2 



(31) 



where C = fin/ a, and finally 

A/" = 



Ti J 

/ damp 



and £ = 



(9 



damp 
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Each line of equation ff28l) can be analyzed to determine what contributions to the resonance 
angle's equation of motion come from the disturbing function's secular and direct portions 
(«Si, iS> 2 , and IZd from equation (12"B"]) ) and what contributions might be expected to arise from 
external damping forces. The first line of equation f[2"5j) contains the immediate contributions 
from the direct part of the disturbing function - that is to say that each term inside the square 
brackets includes a factor ofTZo, but also included are higher order cross term contributions 
to the motion arising from secular effects on Tip. The second line contains the contribution 
from the secular parts of the disturbing function - for the circular restricted three body 
problem, this line vanishes. The final line of equation ( 128]) originates from the presence of 
damping forces like those introduced in equation ([25]) . 



We can expand each operator and determine to what order the eccentricity contributes 
to the overall magnitude in each individual term therein: 



J(aQ) = -C 2 U D 

j(q) = -c 2 n D 



Kd a (a 2 d a Q) - 2a 2 h a h e d e d a Q + a 2 h 2 e d 2 e Q 



+ h a ah e d e Q 



h a ah e d e Q 



q(3* 



a 2 h a d e d a Q - a 2 h e d 2 Q + ( + h a ) ad e Q 



h 2 a a 2 d 2 a Q - 2a 2 h a h e d e d a Q + a 2 h 2 e d 2 Q - 
-aC 2 e- 1 S 2 



JC(Q) = -C 2 e- l S 2 



a 2 h a d e d a Q - a 2 h e d 2 Q +(^r + 0}ad e Q 



(32) 
(33) 

(34) 

(35) 



J(aQ) 



C 



a 



2 



J{Q)=C 



3 

2 ~, 



2 r 



(36)" 



2 ~, 



2 ~, 



Ma 2 h a d 2 a Q+ \ Eeah a - -Afa 2 h e J d e d a Q - Sq(3d 2 e Q + a \SqP* + ^KJ d e Q 

(37) ' 

It is important to consider the magnitude of each of these operators to determine their 
significance in the equation of motion given in equation ( )28|) . Figures [TOj [TTJ and [T2l show 
the resulting magnitudes as a function of the test body eccentricity e for the terms between 
the square brackets in line 1, 2, and 3 of equation (1281) . respectively. 



3.2. Analysis of the Expansion Terms 



In this section we focus on the sin part of equation (128]) - the first line of the equation 
of motion for the resonance angle. As written in this compact form, there are only five 
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Fig. 10. — Magnitudes of coefficients for terms appearing in the first line (sin0 part) of 
equation (1251) as a function of test body eccentricity e. Solid curves represent positive values, 
and dashed curves negative. The magenta curve is the coefficient to 4> (second term in square 
brackets of line 1 of the equation). The cyan curve is the usual pendulum term (first term). 
The red curve is JilZo), the coefficient to cos0 (last term). The blue curve is J(aS2), the 
coefficient to cos (Aw) (fourth term). The black curve is J(aS{) (third term). The figure 
shows that, for low values of eccentricity, the pendulum term is not necessarily the dominant 
term. The figure depicts values for 2:1 external resonances, with n = 2tt, fi = 10~ 3 , and 
e p = 0.1. 
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Fig. 11. — Magnitudes of coefficients for terms appearing in the second line (sinAcc part) 
of equation (1251) as a function of test body eccentricity e. Solid curves represents positive 
values, and dashed curves negative. The magenta curve is the coefficient to w (first term in 
square brackets of line 2 of the equation). The red curve is K,(1Zd), the coefficient to cos0 
(last term). The blue curve is K.(aS2), the coefficient to cos (Aw) (third term). The black 
curve is JC(aSi) (second term). The only term independent of additional contributions from 
the angles (0 and Aw) in this figure is the black curve This figure depicts values for 2:1 
external resonances, with n = 2tt, \i = 10~ 3 , and e p = 0.1. The relative magnitudes are 
dependent on different factors of e p . All strengths vanish for e p = 0, i.e., for the case of the 
circular restricted 3-body problem. 
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e 



Fig. 12. — Magnitudes of coefficients for terms appearing in the damping terms (third line 
of the equation) of equation f l28l) as a function of test body eccentricity e, e.g., damping 
terms that may be considered as contributions to h and e (in addition to the usual Lagrange 
planetary equations) as in equation ( |25l) . Solid curves represent positive values, and dashed 
curves negative. The magenta curve is the coefficient to (the first term). The red curve 
is J(1Ze>), the coefficient to cos0 (fourth term). The blue curve is J7"(a5 2 ), the coefficient 
to cos (Aw) (third term) - it is not present in the circular restricted 3-body problem. The 
black curve is J(aS{) (second term). The cyan curve is Mpn p . This figure depicts values 
for 2:1 external resonances, with n = 2tt, fi = 10~ 3 , and e p = 0.1. Damping parameters are 
r a = 10 4 years and r e = 10 3 years, where we have assumed a simple definition for damping 
parameters a/a = —t^ 1 and e/e— —t~ x . 
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contributing terms contained within the square brackets. The first term, A/"pn p sin </>, when 
considered on its own, has been the focus of intense study in the circular restricted three 
body problem, and is called the pendulum model for obvious reasons. This model can be 
generalized by expanding each term within the square brackets via equations ( 15Tp . ( 11%]) . f[2"T|) . 
( 132|) . and (|33|) and keeping only the lowest eccentricity order contributions from each. Doing 
so gives the result 

~ _ sin 0[ _ 3/2/i Q C(pn p + 0) Ji e 9 - g 2 CJi0 e 9 " 2 + 2C 2 qh a T l {2ad a f sA + f s>1 ) e q 
+ C 2 g 2 e p / S)2 Ji cos (Aw) e 9 " 3 - C 2 q\q - 2) J? cos0 e 2{9 " 2) ] . 

For non-circular perturber orbits, there is an additional contribution proportional to e 2 that 
has been excluded here, which may be significant given large enough values of e p . Notice that, 
depending on the particular value of q, there are up to four different orders of eccentricity 
that show up, so we collect them in eccentricity order, 

OFie" {2qh a C (2ad a f sA + f sA ) - 3(p + q)pn p - (3(p + q) + q 2 e~ 2 ) </>} 
+ C 2 q 2 e q - 3 J r i {e p f S)2 cos (Azj) - q{q - 2)e 9 ~ 1 J"i cos <p[ 



smi 



(39) 



However, the number of terms shown here can be further reduced. Although we are interested 
in a lowest eccentricity order study, we must also consider the lowest order terms in C oc fi 
as well. For a 2:1 period ratio, test body mean motion n = 2tt, mass ratio fi ~ 10~ 3 , and 
perturber eccentricity e p = 0.1, Figure HO] depicts the corresponding relative strengths of 
each term contributing to the sin0 part (first line) of equation ( l28lh 



For eccentricities e ~ 0.1, two terms in the expansion dominate. One term is the 
usual pendulum term, but the other term is something like a modified damping term (i.e., 
a term multiplying the time derivative <fi). Although, |0| is generally less than unity for 
instances of libration, contributions from the term can become quite large for sufficiently 
small eccentricities due to its e _1 dependence. In the figure, contributions from <fr and the 
pendulum term have the same sign (both curves are dashed lines, meaning less than zero 
on the log-log plot). To achieve nodding by way of the <p term, one requires <p < in 
addition to satisfying the constraint on the upper bound of the eccentricity, e < 0.1. The 
condition that <fi < for nodding to occur provides a simple explanation for the observation 
that the resonance angle tends to circulate on average over secular timescales in a preferred 
direction, usually resulting in a graph of <fi versus t which is reminiscent of a stair case. The 
figure actually shows that this term is comparable to the pendulum term for all values of 
eccentricity, however it works to counteract the pendulum term for instances where > 0, 
and has opportunity to overwhelm the pendulum term only for \<p\ > 2. 

Another term of significance appears in equation fl39|) . For sufficiently small test body 
eccentricities, the term that is proportional to e~ 2 cos </> plays an important role. As shown 
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in Figure [101 the cos coefficient (red curve) surpasses the pendulum term (cyan curve) for 
eccentricities lower than ~ 0.03. The greatest lower bound on eccentricity required for this 
term to take effect is 



When evaluated using the typical value \i = fO -3 in the external 2:1 case (J-"i ~ 0.27), the 
above expression gives the restriction e < 0.033. Note that the e -2 cos0 term is not present 
for those resonances with q = 2 (e.g., 3:f , 5:3, 7:5, etc.), and has opposite sign for resonances 
with q = 3 (e.g., 4:f, 5:2, 7:4, etc..) and above. However, for resonances of higher rank 
(higher q values), the eccentricity order increases relative to the order of the pendulum and 
the terms, diminishing its importance. 

Keeping only the two largest terms from Figure [TUl the equation of motion can be 
reduced to the form 



This equation exhibits some of the nodding behaviors we see in the full problem (see sec- 
tion H]). The parameters and initial values must be fine-tuned for cases where eccentricities 
(and periastra) are independent of time (i.e., only special values allow for nodding). Under 
conditions where the eccentricity has time dependence, however, nodding is a robust phe- 
nomenon^ is much easier to find parameters for which nodding occurs). For the sake of 
simplicity, we adopt the parametric description for the time dependence of eccentricity 



In many instances of nodding found from the numerical studies in section 2 but not featured 
in the figures, the test body's orbital eccentricity evolves through large swooping double 
arches, spanning several orders of magnitude, reaching values as high as a few times f0 _1 
and as low as fO -2 — f0~ 3 (see top panel of Figure [6]). This double arched pattern exhibited 
by the test particle's orbital eccentricity is governed by the secular time scale, which typically 
falls in the range ~ fO 2 — fO 3 libration times, but can be much longer (e.g., Michtchenko et 
al. 2008b). The ansatz of equation ( 14 ip is used in order to model this behavior. 

To demonstrate that this model exhibits some nodding behaviors, we integrate equation 
( 140|) using an adaptive fourth order Runge-Kutta integration scheme with e = constant and 
with e(t) given by equation ( )4"Tj) . We note that this is a simple study in which the eccentricity 
evolution being used here is totally prescribed without feedback from the specific orbital 
parameters, so we do not expect the model to exhibit all of the intricate details exhibited 
by the full 3-body simulations presented in section |2j For the sake of definiteness, we take 




~ CJ-\ sin 6n p e + e 1 (f> — CT\t 2 cos 



(40) 



e(t) = e - e a sin 2 (a; e t) . 



(41) 
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n ~ 1, C ~ 10 -2 , and u e ~ 10 -3 yr -1 - the precise values of these parameters will depend on 
the orbital angles and mean motions, and the values we choose corresponds to orbital periods 
on the order of years, not days. Under the parametrically evolving eccentricity of equation 
(|4ip . nodding is a robust phenomenon, even when considering the pendulum term alone. 
Guided by the pendulum model of the circular restricted 3-body problem, we take the initial 
conditions (<j), (ft) = (e, 0), which places the system in a dynamically vulnerable position near 
a separatrix - any perturbations that supply additional action should eventually cause the 
pendulum to circulate rather than oscillate. Figure [T3] demonstrates how a small amount of 
eccentricity variation can send such a system into bouts of nodding. However, the system 
need not be prepared in such a dynamically sensitive manner to see nodding. Figure IT4l 
shows a system that would be stable in the absence of variable eccentricity, but where the 
inclusion of sufficient eccentricity cycling provides enough added action to induce nodding. 
In both Figures [TBI and [T4l we use T\ — 0.27, which is the approximate value obtained from 
equation ([271) corresponding to a 2:1 period ratio with a ~ 0.63 (see MD99). Figure [151 shows 
the results when taking T\ = 1 and using various combinations of the terms in equation (HUl) 
along with a time varying eccentricity as defined in (JUJ. This figure demonstrates that 
the pendulum model alone is not enough to recover nodding behaviors, but nodding does 
appear in models that combine the pendulum term with either of the two additional terms in 
equation (T4T)]) . Taken together, Figures [TBI - [T5l show that nodding behavior arises naturally 
in modified pendulum equations, such as those resulting from the expansion of the previous 
section. 

It may be worth noting that another term in equation (IB9"j) is of order e~ 2 , which was 
ignored due to its contribution from e p . However, this term could be of importance in cases 
where the test particle's eccentricity is very low (e <C e p ). Including this term leaves us with 
a form dependent upon (Aw), 



CJ-j sin . 



6n p e + e — CT\e cos — Cf s ^e p e cos (Aw) . (42) 



This equation will exhibit nodding, even for constant eccentricity (no time variation), as 
long as the eccentricity is sufficiently low and the periastron circulates on secular timescales. 

We caution that this model certainly does not exhaust all possible terms that can lead 
to resonance angle nodding. In the regime where the perturbing planet's eccentricity is non- 
negligible e p > 0.1, terms of first order in the planet's eccentricity can become significant. 
In particular, the term in the disturbing function that goes like ~ e p cos(0 + (Aw)) can 
lead to some interesting modulations in the resonance angle, which ultimately may drive the 
resonance angle's phase into the vicinity of a separatrix and induce nodding (K. Batygin, 
private communication). 
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Fig. 13. — Solutions to the pendulum equation with constant eccentricity (black curve) 
and with eccentricity that varies by 0.1% on a secular timescale (blue and red curves). 
The system is set in a dynamically sensitive state near a separatrix with = 0.001 and 
= and integrated using 4 th order Runge-Kutta scheme for 16 full secular cycles. The 
eccentricity is initially small (e = 0.001), and the pendulum term has no knowledge of the 
inner body's eccentricity. The introduction of time varying eccentricity (even with small 
amplitude) induces circulation in the motion of the pendulum. The eccentricity cycle is 
phase shifted by 7r between the red and blue curves, where the blue curve begins at peak 
eccentricity. The bottom panel shows the phase trajectory for the resonance angle solutions 
in the top panel. The resonance angle is on phase trajectory that is very near a separatrix of 
the phase space. Note that in the phase diagram, the black curve is confined to < < 2n, 
the blue curve is confined to < d> < 6n, and the red curve is confined to — 12n < d> < 107T. 
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Fig. 14. — Solutions to the pendulum equation for initial values e = 0.1, <fi = 1.2 n, and = 0. 
Both solutions presented here are for the pendulum term with (i) constant eccentricity (blue 
curve) and (ii) with eccentricity that varies by 99% on a secular timescale (black curve). 
When eccentricity is constant, the system remains in an oscillatory state with amplitude 
A(f) « 7r/4 around the pointy ~ n. For time varying eccentricity, the librations of the 
resonance angle become unbound before 2 complete secular cycles, and nodding ensues. The 
bottom panel shows the phase trajectory for the resonance angle solutions in the top panel. 
Oscillations of the resonance angle occur around points that are equivalently located at = 7r 
(modulo 27r). 
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The model developed in this section contains one important deficiency. The inner sep- 
aratrix of the external resonance problem is completely missing from this model along with 
the unique characteristics distinguishing the external resonance from the internal one. This 
separatrix originates because of a bifurcation that occurs for the external resonance when 
e ~ 0.1, and is responsible for asymmetric resonances that have been observed in simulations 
for that regime (e.g., Callegari et al. 2004, Lee 2004). As a result, = ir (or a point nearby) 
becomes a hyperbolic fixed point. Given that the test body's eccentricity must be sufficiently 
large before this bifurcation occurs may suggest that more terms of the disturbing function 
are required in the time averaged treatment to recover these dynamics. More work must be 
done in order to elucidate this issue. 



4. Conclusion 

This paper presents an investigation of nodding behavior for planetary systems that 
are near mean motion resonance, with a focus on resonance angles corresponding to the 2:1 
MMR. For systems that experience the nodding phenomenon, the resonance angle librates 
(oscillates) for several cycles, then circulates for one or more cycles, and then resumes its 
oscillatory motion (libration). The process repeats, so that even though the resonance angle 
is primarily in oscillation, the phase of the resonance angle nonetheless accumulates over 
time. In the extreme version of this behavior, the resonance angle can oscillate for one 
cycle, circulate for the next cycle, and then repeat the process; the resonance angle moves 
continually back and forth between the two types of motion, so that the resonance angle has 
an effective period of oscillation that is ~ 2 times longer than the usual period for MMR 
(see Figures 3, 8, and 9). Our numerical exploration of Section 2 shows that the nodding 
phenomenon arises in a wide variety of planetary systems (see Figures 1-9). 

Nodding can be described as complex motion near a separatrix in the phase space of the 
resonance angle. Both internal and external resonances can exhibit nodding, but there exist 
prominent qualitative differences in the nodding signatures between the two configurations. 
The phase space for internal resonances contain one separatrix, whereas the phase space 
for external resonances can contain two distinct separatrices. The qualitative differences 
are mainly due to the existence of asymmetric external resonances, which arise when the 
orbital eccentricity for the outer (smaller) body becomes sufficiently large. Circulation of 
the resonance angle over secular times is common when the planet's orbital eccentricity is 
sufficiently large, i.e., e p > 0.02. Nonetheless, the essential ingredients for nodding behavior 
are the presence of a separatrix and a means to cross it. Note that separatrix crossing can 
occur in sufficiently complex systems (e.g., the driven, inverted pendulum; Acheson 1995) 
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Fig. 15. — The result of numerically integrating equation ( 140]) using parameters T orb = 100 
days, T sec = T or b x 10 2 94 , \i = 10~ 3 , eo = 0.32, e a = 0.28, J~\ = 1, and the orbiting bodies 
in a perfect 2:1 orbital period ratio. The blue curve shows the solution for the pendulum 
term alone; the black curve shows the solution for the pendulum term plus the term that 
goes like 0; the cyan curve includes the pendulum term and the e -2 cos <f> term; the red curve 
shows the solution when including all terms in equation ( 140)) . The only model out of the 
four depicted here that doesn't exhibit nodding is that consisting of the pendulum term 
in isolation. Inclusion of the two terms of lower eccentricity order and in any combination 
results in nodding. The timescale depicted in the top panel is in years. The bottom panel 
shows the phase trajectory for the resonance angle solutions in the top panel. 
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and can arise due to chaos (e.g., Morbidelli & Moons 1993). 

Exoplanets that transit their host stars are expected to exhibit transit timing variations 
when additional perturbing bodies are present, and the variations are greatest when the 
perturbing body is either in or near mean motion resonance with the transiting body. For 
systems where the perturbing body is near MMR, there is a possibility that the resonance 
angle may undergo nodding. In such systems, both the amplitude and the period of the 
TTVs depend on the window of time over which the system is observed (see Figures E] and 
[9]). If the observations are made over a time interval where the resonance angle of the 
system oscillates for many cycles, the TTVs have their usual interpretation. In the limit 
where the resonance angle goes back and forth between oscillating and circulating every 
cycle, the effective frequency of the TTVs is lower (than in the purely oscillatory case) by 
a factor of ~ 2. Nodding behavior that is intermediate between these two cases is also 
possible, and produces TTVs with intermediate frequencies (see the power spectra in the 
bottom panels of Figures [8] and [9]) . Note that the amplitudes of the TTVs vary with the 
mode of nodding/oscillation (for the same masses). The added complexity in the dynamics 
due to nodding can thus introduce corresponding difficulties in interpreting the source of 
transit timing variations. This possible complexity should be kept in mind when searching 
for hidden exoplanets through measurements of transit times. 

Using Lagrange's planetary equations of motion, we have derived a set of generalized 
equations of motion for the resonance angle by including additional terms in the expansion 
in order to account for nodding behavior (see Section [3]). This derivation uses the time- 
averaged disturbing function and initially keeps terms up to second order in eccentricity. As 
expected, the initial expansion includes a large number of terms. We then performed an 
analysis that uses results from the full numerical treatment to determine the relative sizes of 
the various terms in the expansion over the parameter space of interest (see Figures 10, 11, 
and 12). For some parameter values, we found that some of the "higher order" terms (those 
left out of the MD99 derivation) can dominate over those terms commonly used in deriving 
the pendulum model, which provides a standard description for MMR. 

Given the expected magnitudes of the expansion terms, we have constructed a modified 
model for MMR, where the equation of motion includes two additional terms and thus allows 
for more complex dynamics (see equation [40]) . This equation exhibits the nodding behavior 
found for interior resonances of the full problem. However, nodding only occurs for particular 
values of the eccentricity and argument of periastron, where both variables are considered 
as constants. We then generalized this model one step further by allowing the eccentricity 
to change with time (see equation [H]), motivated by the secular cycles of eccentricity 
variation often seen in two planet systems (including the orbits of Jupiter and Saturn in our 
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solar system - see MD99). With this generalization to include eccentricity cycles, nodding 
is a robust phenomenon and occurs for a wide range of the other parameters. Even when 
considering the pendulum term in isolation, the cycling of eccentricity on secular timescales 
can provide sufficient perturbations to induce nodding behaviors in certain dynamically 
vulnerable configurations. 

The model from section [3] contains at least one deficiency. The inner separatrix of the 
external resonance problem is missing from the model, along with the unique characteristics 
distinguishing the external resonance from the internal one. This additional separatrix orig- 
inates for the external resonance for larger outer body eccentricities, when e ~ 0.1, and is 
reminiscent of the pitchfork bifurcation shown in Figure [71 This separatrix is responsible for 
the asymmetric resonances that have been observed for that regime. As a result, some region 
in the neighborhood of <fi = tt contains a hyperbolic fixed point. Given that the eccentricity 
of the test body must be significantly large before this bifurcation occurs, additional terms 
in the expansion of the disturbing function may be required to recover these dynamics. This 
issue is left for future work. 

The results of this work, and the existence of nodding phenomena in general, have two 
principal implications: [1] Planetary systems near MMR display complex and sometimes 
unexpected behavior, so that nodding poses a rich set of dynamical questions for further 
work. It would be interesting to determine the minimal requirements for a dynamical system 
to exhibit nodding, and to explore the relationship between nodding and chaotic motion. [2] 
The main application of this work to observations lies in the interpretation of transit timing 
variations. If an observed system with TTVs experiences nodding, then the inferred mass and 
orbital elements of the unseen perturbing body could vary, depending on the time interval 
of observation (Figures [8] and [9]). A more complete exploration of parameter space, with 
a focus on the TTV signals (analogous to Veras et al. 2011), should be undertaken in the 
future. 

This paper benefited from discussions with many colleagues, especially Konstantin Baty- 
gin, Philip Holmes, Greg Laughlin, Man Hoi Lee, and Matthew Payne. We also thank an 
anonymous referee for useful comments that improved the manuscript. This work was sup- 
ported by NSF grant DMS-0806756 from the Division of Applied Mathematics, a Rackham 
Predoctoral Fellowship from the University of Michigan (JAK), NASA grant NNX11AK87G 
(FCA), and NSF grants DMS-0907949 and DMS- 1207693 (AMB). 



-44- 



REFERENCES 

Acheson, D.J. 1995, Proc R Soc Lond, A 448, 89, 95 

Adams, F. C, & Laughlin, G. 2006, ApJ, 649, 992 

Adams, F. C, Laughlin, G., & Bloch, A. M. 2008, ApJ, 683, 1117 

Agol, E., Steffen, J., Sari, R., & Clarkson, W. 2005, MNRAS, 359, 567 

Batygin, K., & Morbidelli, A. 2011, Celest. Mech., Ill, 219 

Batygin, K., & Morbidelli, A. 2012, A J, submitted 

Beauge, C. 1994, Celest. Mech., 60, 225 

Beauge, C, & Michtchenko, T. A. 2003, MNRAS, 341, 760 

Bodenheimer, P., Laughlin, G., & Lin, D.N.C. 2003, ApJ, 592, 555 

Brouwer, D., & Clemence, G.M. 1961, Methods of Celestial Mechanics (Academic Press) 
Callegari Jr., N., Michtchenko, T., & Ferraz-Mello, S. 2004, Celest. Mech., 89, 201 
Cochran, W. D., Fabrycky, D. C, & Torres, G. et al. 2011, ApJS, 197, 7 



Exoplanets Data Explorer, September 2012, http://exoplanets.org/ 



Fabrycky, D. C, Lissauer, J. J., & Ragozzine, D., et al. 2012. larXiv: 1202.6328^ 2 
Ferraz-Mello, S. 1988, AJ, 96, 400 

Ferraz-Mello, S., Beauge, C, & Michtchenko, T. 2003, Celest. Mech., 87, 99 
Greenberg, R., & Franklin, F. 1975, MNRAS, 173, 1 
Goldreich, P., & Soter, S., 1963, Icarus, 5, 375 

Guckenheimer, J., & Holmes, P. 1983, Nonlinear Oscillations, Dynamical Systems, and Bi- 
furcations of Vector Fields ( Springer- Verlag New York) 

Henrard, J., Lemaitre, A., Milani, A., & Murray, C. D. 1986, Celest. Mech., 38, 175 

Holman, M. J., Fabrycky, D. C, Ragozzine, D. et al. 2010, Science, 330, 51 

Hut, P. 1981, A&A, 99, 126 



-45 - 

Ketchum, J. A., Adams, F. C, & Bloch, A. M. 2011a, ApJ, 726, 53 
Ketchum, J. A., Adams, F. C, & Bloch, A. M. 2011b, ApJ, 741, L2 
Kley, W., Lee, M. H., Murray, N., & Peale, S. J. 2005, A&A, 437, 727 
Lecoanet, D., Adams, F. C, & Bloch, A. M. 2009, ApJ, 692, 659 
Lee, M.-H., & Peale, S. J. 2002, ApJ, 567, 596 
Lee, M.-H., 2004, ApJ, 611, 517 

Michtchenko, T., Beauge, C, & Ferraz-Mello, S. 2008a, MNRAS, 387, 747 
Michtchenko, T., Beauge, C, & Ferraz-Mello, S. 2008b, MNRAS, 391, 215 
Morbidelli, A., & Moons, M. 1993, Icarus, 102, 316 

Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge 
Univ. Press) (MD99) 

Nesvorny, D., & Morbidelli, A. 2008, Icarus, 688, 636 

Quillen, A. C. 2006, MNRAS, 365, 1367 

Rein, H., & Papaloizou, J.C.P. 2009, A&A, 497, 595 

Snellgrove, M. D., Papaloizou, J.C.B., Nelson, R. P. 2001, A&A, 374, 1029 
Veras, D., Ford, E. B., & Payne, M. J. 2011, ApJ, 727, 74 



This preprint was prepared with the A AS IATgX macros v5.2. 



